Pairing in a two-dimensional Fermi gas with population imbalance 
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Pairing in a population imbalanced Fermi system in a two-dimensional optical lattice is studied 
using Determinant Quantum Monte Carlo (DQMC) simulations and mean-field calculations. The 
approximation-free numerical results show a wide range of stability of the Fulde-Ferrell-Larkin- 
Ovshinnikov (FFLO) phase. Contrary to claims of fragility with increased dimensionality we find 
that this phase is stable across wide range of values for the polarization, temperature and interaction 
strength. Both homogeneous and harmonically trapped systems display pairing with finite center 
of mass momentum, with clear signatures either in momentum space or real space, which could be 
observed in cold atomic gases loaded in an optical lattice. We also use the harmonic level basis in the 
confined system and find that pairs can form between particles occupying different levels which can 
be seen as the analog of the finite center of mass momentum pairing in the translationally invariant 
case. Finally, we perform mean field calculations for the uniform and confined systems and show the 
results to be in good agreement with QMC. This leads to a simple picture of the different pairing 
mechanisms, depending on the filling and confining potential. 
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I. INTRODUCTION 

The question of pairing in polarized fermionic systems 
came to the fore shortly after superconductivity in un- 
polarized systems was explained by BCS pQ as being 
due to the formation of Cooper pairs with zero center 
of mass momentum. Fulde and Ferrell [2] and inde- 
pendently Larkin and Ovchinnikov [3] proposed similar 
but not identical mechanisms whereby the fermions form 
pairs with nonzero center of mass momentum. We will 
refer to such a phase as the FFLO phase. On the other 
hand, Sarma [3] proposed a mechanism where, in spite 
of the mismatch in the Fermi momenta due to the spin 
population imbalance, the pairing occurs with zero center 
of mass momentum. Verifying these predictions experi- 
mentally proved difficult in condensed matter systems [5] . 
However, thanks to rapid experimental progress in the 
domain of ultra-cold atoms, it is now possible to study 
such population imbalanced systems. Fermionic atoms 
are made to occupy two hyperfine states thus emulating 
a system with "up" and "down" spins. An advantage of 
these systems is that the population imbalance (the po- 
larization) and the interaction strengths are highly tun- 
able. Such experiments have been performed in three- 
dimensional [HI E] and one-dimensional [5] systems. 

It is by now widely accepted that at T = the FFLO 
phase is robust over a wide range of parameters in one- 
dimensional systems with imbalanced fermion popula- 
tions. This was shown in various numerical studies using, 
for example, QMC 0HU] and DMRG QIHI3]. In a pre- 
vious work we also showed that the FFLO phase is sta- 
ble over a wide range of parameters in the temperature- 
polarization (TP) phase diagram [IB]. This exotic pair- 



ing occurs both in homogeneous and confined systems, 
and has been shown to survive up to relatively high tem- 
peratures (T/Tp sa 0.1) which are achievable in current 
experiments. 

The question of the stability of this phase in higher di- 
mensions remains a subject of debate. It is believed that 
"nesting" of the Fermi surfaces stabilizes FFLO pairing. 
For example in one dimension one wave- vector connects 
all points on the Fermi surfaces of each species, which 
would enable all particles from the Fermi surfaces to par- 
ticipate in the formation of pairs with finite-momentum. 
The effect of "nesting" is considerably weaker in higher 
dimensions. In a two-dimensional lattice system, the 
shape of a Fermi surface depends on the filling. At half 
filling the Fermi surface becomes a square and touches 
the edge of the first Brillouin zone (Van Hove singular- 
ity) . Around this filling, matching of the Fermi surfaces 
becomes more efficient, in other words the "nesting" is 
enhanced as compared to the situation when both Fermi 
surfaces are circular (low filling) . This reasoning leads us 
to expect that FFLO pairing should be more prevalent 
around half filling than at lower fillings. This lattice en- 
hanced stability of FFLO was studied using mean field 
(MF) methods in [19] and [20]. In the latter the authors 
point also at Hartree corrections and domain wall forma- 
tion as additional reasons for enhancement. 

Numerous theoretical studies of the system in higher 
dimensions do not offer a clear conclusion on the stability 
of the FFLO mechanism. In a variational MF study of 
a three dimensional system in the continuum with and 
without a trapping potential it is observed that FFLO 
is a fragile state which can be realized only in a tiny 
sliver of the interaction-polarization phase diagram |21j . 
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Furthermore, this study showed that in a trap, FFLO can 
exist only in a thin shell of the atomic cloud. Another 
study of a three dimensional Fermi gas at unitarity [2"2"] 
shows that this phase is competitive over a large region 
in the phase diagram. However, the trap would need 
to be adjusted to allow FFLO to occupy a large enough 
spatial region to be observed. On the other hand, in a 
Bogoliubov-de Gennes study [T3] of a trapped system, 
the calculated phase diagram indicates that the ground 
state of the system is always FFLO for any imbalance up 
to some critical value. 

The unsettled status of this phase in higher dimension- 
ality may be clarified with exact numerical simulations. 
However, simulations of the Hubbard model in three di- 
mensions are not feasible for large systems at low enough 
temperatures due to the severity of the "fermion sign 
problem". On the other hand, exact QMC simulations 
in two dimensions are feasible but so far none have been 
done. In addition, two dimensional systems are interme- 
diate between one dimension where MF is almost certain 
to fail and three dimensions where MF is more reliable. 
Consequently, there has been a concerted, yet incon- 
clusive, effort to understand FFLO physics theoretically 
in two-dimensional systems. Homogeneous and trapped 
two-dimensional polarized Fermi gases have been stud- 
ied with MF calculations which exclude the possibility of 
FFLO pairing (for e.g. [24] and [25]). An interaction- 
polarization phase diagram is shown in |27] where FFLO 
pairing is seen to occupy a wide region. Koponen et al. 
|19j obtain MF phase diagrams in the polarization versus 
filling plane for one-, two- and three-dimensional systems. 
In the two-dimensional system there is a very strong fea- 
ture around the van Hove singularity of the majority 
component and the FFLO phase is present over a wide 
range of parameters around this value. They also show 
temperature-polarization phase diagrams of one dimen- 
sional system which were shown not to agree with exact 
QMC results [16]. The temperature-polarization phase 
diagram in three dimensions is shown as well but not the 
two dimensional case. Studies of quasi two-dimensional 
systems have been done using MF and they predict a 
first-order transition to FFLO at finite temperature [26] . 
Another mean field study of two-dimensional two-orbital 
Hubbard model with p-orbitals and highly unidirectional 
hopping shows enhancement of the FFLO region in the 
phase diagram due to the one-dimensional character of 
the Fermi surface [33J . A DMRG study of population im- 
balanced Fermi gas on two-leg ladders has found FFLO 
pair correlations [2"5] , 

In this paper we present a Determinant QMC (DQMC) 
|29j study of the two dimensional Hubbard model with 
imbalanced populations of up and down spins. In sec- 
tion II we present the model and discuss our results 
for the uniform system in section III. Our main result 
here is the demonstration of the robustness of the FFLO 
phase and the determination of the phase diagram in the 
temperature-polarization plane at low filling. We also 
compare the behavior of the system at low and half fill- 



ings. In section IV we examine the system in a harmonic 
trap. Our conclusions are in section V. 



II. MODEL AND METHODS 

The system of interest is governed by the two- 
dimensional fermionic Hubbard Hamiltonian 

<i,j>a i 




+^T^(j-j c ) 2 (n i i+n j2 ) (1) 
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Where c\ a (c- ) create (annihilate) a fermion of spin 

a = 1, 2 on lattice site i and hi a = c t a c i a is the corre- 
sponding number operator. The near neighbor, < i,j >, 
hopping parameter is t which we take equal to unity to 
set the energy scale. We consider only on-site interaction 
with an attractive coupling constant U < 0. The number 
of particles in each population is governed by its chemi- 
cal potential (/i<t). The harmonic trap is introduced via 
the Vt term in the Hamiltonian where j c is the position 
of the center of the trap (also middle of the lattice) . All 
simulations are performed with periodic boundary con- 
ditions. In the confined case we ensured that the density 
vanishes at the edge of the lattice. 

The main quantities of interest in this study are the 
single particle Green functions G a and the pair Green 
function, G pa i r , 

G*(l) = <ct +iff O, (2) 
G pair (0 = (At +J A;), (3) 
A, = c i2 cn, (4) 

where A i creates a pair on site i. The Fourier transform 
of G a gives the momentum distribution of the spins-a 
species while the transform of G pa ; r yields the pair mo- 
mentum distribution. In the trapped case, the density 
profiles of the two species are also studied. 

We studied this system numerically using the DQMC 
[2"5] algorithm. In this approach, the Hubbard- 
Stratonovich (HS) transformation is employed to decou- 
ple the quartic interaction term into two quadratic terms 
coupling the number operator of each species, n a (i), to 
the HS field which effectively acts as a site and imaginary 
time dependent chemical potential. The fermion opera- 
tors can now be traced out leading to a partition function 
in the form of a product of two determinants, one for each 
spin, summed over all configurations of the HS field. For 
U < and equal populations, /ii = jj,2, the determinants 
are identical; their product is always positive. But in the 
imbalanced case, /ii ^2, the two determinants are no 
longer equal and their product can, and does, become 
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negative leading to the known "fermion sign problem". 
This is the main obstacle to the simulation of this sys- 
tem. We found that at low total filling the sign problem 
is manageable even at large polarizations and low tem- 
peratures. This was not the case closer to half filling. 
Typical simulations of the harmonically confined system 
at low temperature took about two weeks on a 3 GHz 
processor. 

In the presence of the trapping potential, we have also 
studied the system using a mean-field approach. Starting 
from the full Fermi-Hubbard Hamiltonian ([I]), one can 
derive the mean-field Hamiltonian: 



h mf = ipiMi/j + ^Y1 A * A * - E 



' ~ 1 -A? -h 



U 

A,; 
ji2 



(5) 



where A* 
*t = (.. 



U(c\ lC . 



are on-site pairing amplitudes; 
is the Nambu spinor. The 



matrix h depicts the one particle Hamiltonian, namely 



hopping terms between nearest neighbors hi 



-t i £ 



j and chemical potential terms hjj a — ~Pj a = —l^a + 
Vr (j — j c ) 2 fij a- To account properly for spatial inho- 
mogeneities, the BCS order parameter at each site, Aj, 
is an independent variable [3H1 [32 133]: whose value is 
determined, for a given temperature, by a global min- 
imization of the free energy, F = — -g-ln(Z) associated 
with the mean-field Hamiltonian: 



F=--Vln(l 



^A*A 4 -5> 2 , (6) 



U 



where the Afc are the 2N eigenvalues of the Nambu ma- 
trix M; N is the number of sites. The minimization of 
the free energy is performed using a mixed quasi-Newton 
and conjugate gradient method; additional checks were 
performed to ensure that the global minimum has been 
reached. 



III. HOMOGENEOUS SYSTEM 

To set the stage, we start with the homogeneous two 
dimensional Hubbard model with balanced populations 
at total density p = p±+p 2 = 0.3. With balanced popula- 
tions, the pairs form with zero center of mass momentum 
and a sharp peak in the pair momentum distribution is 
expected at k — 0. Figure [l] shows the momentum distri- 
butions for a system with U = — 3.5t, p\ + p 2 = 0.3 and 
f3 = 30 in a 16 x 16 optical lattice. The single particle 
momentum distribution, identical for the two spins, is 
shown in Fig. [l|a) while (b) shows the pair momentum 
distribution. As expected for weak to moderate values of 
\U\, the single particle distribution has the usual Fermi 
form and the pair momentum distribution exhibits a very 
sharp peak at k = indicating pairing with zero center 
of mass momentum. 




FIG. 1: (a) Single particle momentum distribution, 
ni(k x ,k y ) (the same as U2(k x , ky)). (b) Pair momentum dis- 
tribution, n pa i r (k x ,k y ) exhibiting a sharp peak at zero mo- 
mentum. The total density is pi+p2 = 0.3 (pi = p2), ft = 30, 
U — —3.5t and the system size is 16 x 16. 



We now examine the polarized system. To this end, 
the chemical potentials p\ and pi are made unequal so 
that pi ^ pi but p = pi + p2 remains constant. This re- 
quires tuning the chemical potentials appropriately. The 
polarization, P, is defined by 



P : 



Nt -N 2 
N 1 +N 2 " 



(7) 



where N% and N 2 are the total populations of the two 
species. 

Figure [2] shows the momentum distributions for a sys- 
tem with U = -3.5t, P = 0.6, p = 0.3 and /3 = 10 
in an optical lattice of size 16 x 16 for (a), (b) and (c) 
and 10 x 30 for (d). Panels (a) and (b) show the mi- 
nority and majority single particle momentum distribu- 
tions, n\(k x , k y ) and n2(k x ,k y ), respectively. They ex- 
hibit usual Fermi-like distributions. However, the pair 
momentum distribution, n va \ r (k Xl k y ), is strikingly differ- 
ent from the balanced case: It has a volcano-like shape 
with the maximum of the distribution at the rim of the 
crater of radius |fc| = \kp 2 — kp\\. kp\ and kp 2 are the 
minority and majority Fermi momenta respectively. 

In two dimensions, the Fermi surface geometry changes 
with the filling. The behavior exhibited in Fig.[2]is for low 
filling where the Fermi distributions of both species have 
cylindrical shape and the pairs are formed with equal 
probability in all radial directions. In this density regime, 
the signature for the FFLO phase is the presence of a cir- 
cular ridge in the pair momentum distribution as seen in 
Fig. J2jc) . Studying the system in the low density regime 
is interesting because it also approximates the continuum 
conditions. 
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FIG. 2: Momentum distributions of (a) minority and (b) 
majority populations, (c) shows the pair momentum distri- 
bution. The parameters are p = p 1 + P2 = 0.3, P = 0.6, 
j3 = 10, U = —3.5t in an optical lattice of size 16 x 16. (d) 
The pair momentum distribution for the same system but for 
a lattice of size 10 x 30. 



To study possible finite size effects, we performed our 
simulations for systems of various sizes. In particular, 
Fig. [2^d) shows the pair momentum distribution for the 
same parameters as (a,b,c) but with a system of size 
10 x 30. It is seen that the peak in the pair momentum 
distribution is at the same values of \k\ = \kpi — kp2\ as 
the 16 x 16 system. 

We now examine the effect of temperature on the 
FFLO phase. In particular, we map out the phase di- 
agram in the temperature-polarization plane. Thermal 
effects are very important in experiments due to the dif- 
ficulty in cooling fermionic atoms. The inset in Fig. [3] 
shows two-dimensional cuts in the three-dimensional pair 
momentum distribution for a 16 x 16 system with U = 



— 3.5t, p — 0.3 and P = 0.55. We see that as the tem- 
perature is increased ((3 decreased) the FFLO peak at 
nonzero momentum decreases and, in fact, shifts towards 
zero momentum. Our criterion for the appearance of the 
FFLO phase is when the peak of the pair momentum 
distribution is no longer at zero momentum. The ques- 
tion is then what replaces the FFLO phase: have the 
pairs been broken by thermal fluctuations or has the sys- 
tem been homogenized, resulting in a uniform mixture 
of pairs and excess unpaired particles of the majority 
population? The double occupancy, D = (ni{r)ri2(f)) , 
offers a measure of how tightly bound the pairs are: In 
the absence of pairing, D — p\pi while when the pairing 
is complete, D = pi where p\ is the minority popula- 
tion. These limits suggest the use of a normalized form, 
(D— P1P2) / {pi— P1P2), which is now bounded by and 1. 
Note that p\ — Ni/L 2 while (ni(f)) is the average num- 
ber of type 1 particles at r. In the absence of pairing the 
two quantities coincide. We see in Fig. [3] that for f3 > 3 
the normalized double occupancy is essentially constant 
signaling the continued presence of pairs. This means 
that when the FFLO peak first disappears at 4 < (3 < 3, 
the pairs are still formed. We conclude therefore that 
the system leaves the FFLO phase to enter a polarized 
paired phase (PPP) phase. 

When the thermal energy, T = 1/(3 is of the order 
of the pair binding energy, \U\, the pairs are expected 
to break. We see in Fig. [3] that the double occupancy 
decreases precipitously only for (3 < 1 which is consistent 
with the value of 1/|(7| = 1/3.5 in our simulation. Similar 
behavior was found for the one-dimensional system |16j . 
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16x16, P=0.55, p=0.3, U=-3.5 
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FIG. 3: (Color Online). The normalized double occupancy as 
a function of inverse temperature, ft for p = 0.3, P = 0.55 and 
U — —3.5t. The lattice size is 16 x 16. Inset: Behavior of the 
pair momentum distribution as the temperature is increased 
(/? is decreased). 

Note in Fig. [3] that the double occupancy increases 
just before it drops signaling the breaking of the pairs. 
This increase can be understood physically as follows. 
As the temperature is increased, the Fermi distribution 
near the Fermi momentum gets rounded but for \k\ < 
\kp\ the distribution remains saturated. This means that 
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pairing can happen only near the Fermi surface while 
inside the Fermi sea the particles are still blocked by 
the Pauli exclusion principle. Eventually, as T continues 
to increase, the occupation of momentum states inside 
the Fermi sea drops, rather suddenly as shown by our 
simulations, which makes available for pairing a larger 
number of particles causing the double occupancy to rise. 



1 

0.8 
0.6 
0.4 
0.2 
0, 




3.5 - 



>p=0.3, lattice size 16x16,T F =1. 
. p=0.3, lattice size 20x20 
i half filling, p=1 .0, lattice size 1 6x1 6 







0.1 



0.2 0.3 0.4 0.5 



0.6 0.7 



FIG. 4: (Color Online). Finite temperature phase diagram of 
the system at p — 0.3 

The phase diagram is mapped by fixing the polar- 
ization, P, and increasing T until the peak in the pair 
momentum distribution shifts to zero momentum (inset 
Fig.|3|. The phase diagram for p = pi + p2 — 0.3 (circles) 
and p — 1 (squares) is shown in Fig. |4j The solid circles 
show the boundary of the FFLO phase; the open circles 
indicate the largest P at which we were able to study 
the system. Up to these high polarizations the system 
remained in the FFLO phase. The FFLO phase bound- 
ary at low P appears to extrapolate to P ^ as T — > 
for the 16 x 16 system. However, this is an effect of the 
coarseness of the lattice grid. As P decreases, the peak 
in the pair momentum distribution falls between and 
2t:/L and gives the impression of peaking at zero momen- 
tum. The solid triangles show the phase boundary for a 
20 x 20 system; we see that effect is corrected for a while, 
but then even larger systems are needed. This is not pos- 
sible because as T decreases the sign problem becomes 
too severe. We believe that as soon as the system is po- 
larized it goes into the FFLO phase if T is low enough. 
The long dashed line connecting this FFLO boundary 
to the origin schematizes this. Outside the FFLO phase 
the system is in the polarized paired phase (PPP) since 
the pairs are still formed and break only at higher T 
than shown in the figure. The squares in Fig. [4] show 
the phase boundary at these temperatures for the case of 
p = 1 (discussed below). 

It is important to emphasize here that, in our discus- 
sion, the FFLO state is characterized by the behavior of 
the pair momentum distribution: If the peak is at non- 
zero momentum the system is in the FFLO phase. The 
question naturally arises as to whether the FFLO pairs 
have phase coherence and are, consequently, superfluid. 



In the balanced case, the phase diagram in the temper- 
ature versus filling plane was determined for U = — 4 in 
Ref. [IB] • By studying the pairing susceptibility as a func- 
tion of T as in Ref. [15] , we find that in the balanced case 
of our system with U = —3.5, the critical temperature 
is T c s» 0.1 in good agreement with the U = —4 results 
[18] . However, studying the same pairing susceptibility 
in the polarized case showed no sign of s-wave superflu- 
idity in the temperature range attainable by QMC. Our 
numerical results suppport approximate analytic results 
which indicate that polarization may suppress superfluid- 
ity in the FFLO phase [17]. It is, therefore, currently not 
clear if when T is reduced even further, the FFLO phase 
will become superfluid. We note, however, that the cur- 
rent focus of most experimental measures of FFLO is the 
same non-zero momentum peak on which our simulations 
concentrate. 

The phase diagram, Fig. [4] resembles the one found 
in one dimension |16j and shows that FFLO is very ro- 
bust. The Fermi temperature is calculated as usual by 
considering a balanced ideal system and gives for p = 0.3 
a value T> = 1.881 The FFLO phase at high P survives 
up to T — 0.2Tp while in one dimension [16] at p = 0.25, 
FFLO survives up to T = 0.8T F at high P. So, while 
FFLO is still robust in two dimensions, it is more easily 
destroyed by finite T. This is important to keep in mind 
in experiments. 

In a two-dimensional lattice, the Fermi surface geome- 
try evolves with the filling from closed, rotationally sym- 
metric surfaces for low filling to a square at half filling 
to open surfaces for higher filling. Consequently, pairing 
at finite momentum occurs with different symmetries de- 
pending on the filling. The pairs form with equal prob- 
ability in all radial directions in the case of low filling 
while they form in preferred directions when the Fermi 
surfaces are anisotropic. 

As discussed in the introduction, there are claims that 
around the Van Hove singularity the FFLO pairing could 
be enhanced due to increased nesting. Indeed we observe 
that FFLO is stable over a wider range of temperatures 
and polarizations for p = 1. The squares in Fig. [4] show 
the FFLO-PPP boundary in the half filled case. It is seen 
that the FFLO phase persists to higher T than the low 
density case. However, when compared to Tp — 6.28f, 
FFLO is destroyed for T ~ 0.08Tp as compared with 
T m 0.22V f° r the half filled case in one dimension. 

When the populations are imbalanced around half fill- 
ing of the lattice, one can readily see the effect of the 
interaction on the Fermi surfaces in Fig. [5] depicting 
the difference between the Fermi distributions of the 
species calculated using both mean field and QMC meth- 
ods: they almost look like nested squares parallel to each 
other in most of the momentum states, whereas the non- 
interacting ones would look more rounded and not as 
parallel. Similar Fermi surface geometry in the context 
of LO states in 3D have been shown in [30] . 

The reason why the system exhibits such Fermi sur- 
faces can be understood as follows. If we look at the re- 
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FIG. 5: (Color online) Top row: dilference in the momentum 
distributions of majority and minority, showing parallel Fermi 
surfaces from the Mean-Field method (a) and from QMC (b). 
Bottom row: pairing schematic for balanced (c) and imbal- 
anced (d) . In the situation when the populations of fermionic 
species are imbalanced (diagram on the right) a particle from 
the majority species forms a pair with a particle from the mi- 
nority whose Fermi momentum either matches the k x or k y 
coordinate of the majority particle Fermi vector. The pair 
formed has a finite momentum equal to the distance of the 
two Fermi surfaces either along k x or k y . 




gion of k x > we can parametrize the linear part of the 
majority Fermi line as k^ 2 (k x ) = —k x + a-2 for positive 
values and kj 2 (k x ) = k x — cxi for negative values and do- 
ing the same for the minority we have k~^ \{k x ) — ~k x +a\ 
and kj 1 (fc 2; ) = k x — a\. (see Fig. jHJ). Pairing happens 
here for a given k x between the upper part of the ma- 
jority branch and the lower part of the minority branch 
(fey 2 (fc x ) pairs with fe7i(— k x )). The momentum of the 
pair along y is the sum of these momenta and is equal to 
q y = k~j 2 (k x ) + kj^—kx) — a.2 — ol\- Therefore, thanks 
to the parallel Fermi lines, the pairing momentum is in- 
dependent of k x , leading to a strong enhancement of the 
pairing efficiency The same construction can be done in 
the k x direction, matching the y coordinate of the mo- 
mentum vectors and the pairs will be moving along x 
with ±q x . In other words, for each k x , we have, along 
k v , the usual imbalanced ID situation, i.e. two rectan- 
gular Fermi distributions, with different Fermi momenta. 
Again, the crucial point is both the majority and mi- 
nority effective ID Fermi momentum values change the 
same way with k x : the two Fermi surfaces remain always 
at the same distance from each other. This pairing mech- 
anism is illustrated in Fig. [5ji. The excess fermions cor- 
respond to the part of the majority Fermi surface which 
can't be paired this way, i.e. the four regions around 



FIG. 6: Momentum distributions of (a) minority, (b) major- 
ity and and (c) order parameter in k-space calculated using 
the Mean-Field method, p = pi + p 2 = 1. Here P = 0.32, 
P — 25, U = — 3.5t and the lattice size is 79 x 79. The pair- 
ing peaks are symmetric along k x or k y depending on the 
realization. 



(k x = 0,ky = ±7r) and (k x = ±tt, ky = 0). Note that 
in the balanced case, this corresponds to the usual BCS 
pairing on a lattice: a particle of one species from the 
Fermi surface can form a pair with a particle from the 
other species with the Fermi vector of equal length but 
opposite direction (as shown in Fig. |5t). The resulting 
pair has, as expected, a zero center-of-mass momentum. 
The pairing along k x and k y might not seem the most 
intuitive scenario, since one can imagine the pairs form- 
ing with momentum along the diagonal with smaller \k p \. 
Since this pairing was not observed in any of our simula- 
tions, this probably means that, in a mean field approach, 
it only corresponds to a local minimum of the free energy. 
However, since the shape of Fermi surfaces is affected by 
the nature of the pairing, one can not directly compare 
both situations from the present results and a more de- 
tained study is needed, which is beyond the scope of this 
paper. On the contrary, the mean-field simulations show 
sharp peaks either along k x or k y depending on the re- 
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FIG. 7: Momentum distributions of (a) minority, (b) major- 
ity and (c) pairs at p = pi + P2 = 1, obtained from QMC. 
Here P = 0.38, /3 = 10, U = — 3.5i and the lattice size is 
16 x 16. The pair momentum distribution depicts four peaks 
along the k x and k y axis. 



turn, where the momentum of the pair is the difference of 
the Fermi momenta of each involved fermion. When we 
turn to study a harmonically confined system at low fill- 
ing, for which only few harmonic levels are actually filled, 
the translationally invariant momentum space descrip- 
tion is no longer the obvious one. An ideal gas confined 
in a harmonic trap is known to be fully characterized by 
the basis formed by harmonic oscillator wave functions. 
In addition, for low fillings of the lattice, only the bottom 
of the band structure will be filled. Then the kinetic part 
of the Hamiltonian is well described by the free particle 
one with an effective mass m* given by m* = l/2a 2 i, 
where a is the lattice spacing and t the tunneling ampli- 
tude. In the present case, setting the units t = 1 and 
a = 1, the effective mass is therefore m* = 1/2. 

In this section of the paper we explore the description 
of the interacting system in the harmonic basis. This 
transformation is the analog of the Fourier transforma- 
tion used to go from real space to momentum space in the 
case of the free system. We will show that both BCS and 
FFLO models can be translated into the harmonic level 
basis as pairing of particles between harmonic levels and 
look into the limitations of this description. Since we are 
studying a two-dimensional system we use the eigenstates 
of the two-dimensional harmonic oscillator (see for exam- 
ple [38]). Due to rotational symmetry, the n th harmonic 
level is n + 1 times degenerate. We will use the labeling of 
the states as follows: n is the principal quantum number 
and m = —n, —n + 2, n is the orbital angular momen- 
tum quantum number. For simplicity we will sometimes 
use k to label the set of quantum numbers, k — (n,m). 
Taking the normalized Harmonic oscillator wave function 
for a particular level to be $ n ,m(*) (where i is the lattice 
site) we define a creation operator of a particle in a level 
as: 



alization (see Fig. |6| and in the Quantum-Monte Carlo 
simulations, since we average over all realizations, we see 
that the pair momentum distribution exhibits four peaks: 
two along k x and two along k y (see Fig. [7| . It is impor- 
tant to notice a very good agreement between the results 
obtained by MF and QMC methods. Finally, we have 
also observed, as expected, that the value of the position 
of the peaks, i. e. the center of mass of the pairs, increases 
with large population imbalance. 



1 



(8) 



which, in the continuum limit, leads to properly anti- 
commutating fermionic operators. 

We calculate the single particle Green function be- 
tween levels for each species as follows: 



G<t(k, «') 



- /\T/t 



(9) 



IV. HARMONICALLY CONFINED SYSTEM 

One is used to describing free fermions on a lattice 
using intuition built on the free electron model. Each 
particle occupies a state with particular momentum k 
and at T=0 the filled state with the highest k is called 
the Fermi level. BCS pairing mechanism is understood 
as pairing between fermions from the Fermi surface with 
opposite spins and opposite momenta. In this descrip- 
tion the FFLO pairing model predicts forming a pair of 
fermions from different spin species with a finite momen- 



As pairing is our main interest of investigations we also 
define a pair Green function using the creation and an- 
nihilation operators of a pair of fermions. Similarly to 
the homogenous case where the pairs are formed between 
particles having different momenta, the pairs here can 
have constituents occupying different harmonic levels: 



G palr (K, K') = (^1.^1,2^ K ^ K>A ) 



(10) 



In this section we present results for these correlation 
functions obtained using both QMC and MF method of 
balanced and polarized systems with low filling of the lat- 
tice. All QMC results were done on a 20 x 20 lattice at 
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FIG. 8: (Color online). Single particle Green function in the 
harmonic level basis using QMC, in the balanced case. The 
total number of particles is 22.3, i.e w 11 particles per spin. 
As one can see, the single particle Green function value on 
the diagonal sharply drops just before the 5 th level (n = 4) 
corresponding to 10 harmonic states, roughly the number of 
particles per spin. The off-diagonal elements are small com- 
pared to the diagonal ones, emphasizing the accuracy of the 
harmonic description of the system. States are labelled with 
K = (n, m) and only the principal quantum number n is dis- 
played on the x and y axis. 



the inverse temperature /3 = 10 with interaction strength 
U = —3.5 and the trap potential Vt = 0.065 which trans- 
lates to an effective harmonic frequency uj = 0.5. The 
simulations done using the Mean-Field method were per- 
formed on a bigger lattice of 41 by 41 sites, at the in- 
verse temperature of /3 = 25 and taking the interaction 
strength to be U = —3.0 In the figures only the n val- 
ues are explicitly written, but correlations are calculated 
between all different n and m values. The m levels are 
arranged from m=-n to n from left to right (or bottom 
to top). In the balanced case shown in Fig. [8] the sin- 
gle particle green function is mainly diagonal which indi- 
cates that in this regime the harmonic level basis offers a 
good description of the system. The diagonal part is the 
occupation of levels and where it drops to zero one can 
define the Fermi level. We compare these results to those 
obtained using the Mean-Field method. Both single par- 
ticle as well as pair Green function shown in Fig. [T0| agree 
qualitatively to the QMC results. The small off diago- 
nal values in QMC, which are not present in the MF re- 
sults, stem from the exact treatment of the interactions 
in QMC, not taken into account in the MF calculations. 



FIG. 9: (Color online). Pair Green function in the harmonic 
level basis using QMC. The total number of particles is 22.3 
and the populations are balanced. One can clearly see that 
the pairing is maximum at the Fermi level n = 4 — 5, with op- 
posite magnetic quantum numbers m. Off-diagonal pairing, 
e.g. between k — (6,m) and «' = (4, to'), is almost negligi- 
ble. By diagonal pairing we mean pairing between levels with 
equal principal quantum numbers. 



In the regime of much higher fillings of the lattice (for 
example around half- filling) , the effective mass approach 
is no longer valid and the MF results show that the har- 
monic basis is no longer a relevant one. We do not have 
any QMC results in that regime due to the sign problem. 

In the balanced population case, both Fig. |9j (QMC) 
and Fig 10 (MF), emphasize that the pairing is maximum 



around the Fermi level and happens between particles 
from levels with the same n and for opposite m and m! 
values such that the total orbital angular momentum of 
the pair is 0. This situation is similar to the free particle 
case, where the pairing occurs mostly between the +kp 
and — kp states. 

At low imbalance, one observes that the pairing mostly 
occurs between the same levels, for instance in Fig. [TT] 
where one observes diagonal pairing for n — 3 and for 
n = A. However, one observes an off diagonal feature 
appearing that corresponds to pairing between the lev- 
els n = 3 and n = 4. When the system is imbalanced 
even more, the off-diagonal feature becomes the main 
pairing amplitude. For instance, as shown in Fig. |12[ 
corresponding to a polarization P = 0.22, the diagonal 
pairing has almost completely disappeared and the pair- 
ing mostly occurs between the levels n — 3 and n = A. 
Since it corresponds to a pairing between an odd and 
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FIG. 10: (Color online). Single particle and pair Green func- 
tion in the harmonic level basis using MF. Total number of 
particles is 80.4 and the populations are balanced. As in 
FigjSJ the single particle Green's function is diagonal, with 
a value equal to 1 up to the Fermi level (n « 8) dropping to 
zero after. The pair Green's function emphasizes the diagonal 
pairing (n, m) o (n, — m). 



FIG. 11: (Color online). Single particle (a) and (b) and pair 
Green functions (c) in the harmonic level basis (QMC results) 
for a low polarization situation (P=0.11). The total number 
of particles is 25.5. Even though the Fermi-levels between the 
two species no longer match, the pairing is still diagonal for 
n — 3 and for n = 4 levels. However, one observes an off 
diagonal feature appearing that corresponds to pairing be- 
tween the levels k = (4, —4) and «' = (3, 3) and respectively 
k — (4, 4) and k! — (3, —3) as indicated by arrows. 



even level, it is impossible to match the m values and 
get the pairing with total angular momentum zero. We 
observed that the strongest pairing happens, for exam- 
ple, between re — (4, —4) and re' — (3, 3) and analogously 
between re = (4,4) and re' = (3,-3). There is, in addi- 
tion, a small contribution from the levels re = (4, —2) and 
re' = (3, 1) and re = (4, 2) and re' = (3, —1). In both cases 
the sum of the orbital angular momentum is non-zero. 
Imbalancing the system even more we arrive at the sit- 
uation where the difference between the Fermi levels of 
each species is np2 — npi =2. As illustrated in Fig. [13] for 
P=0.37 the pairing occurs between the levels n — 5 and 
n = 3 and also n — 4 and n = 2 which means that the 
system can now achieve pairing with zero total orbital 
angular momentum. Still, there is small contribution of 
pairing between re = (5, —5) and re' = (3, 3) and re = (5, 5) 
and re' = (3,-3), for which Am = ±2. For a compari- 
son we show the results from the mean-field simulations, 
depicting a similar behavior. In the realization shown in 
Fig. [14] the Fermi levels of each species are n = 7 and 
n = 9 and we can see the pairing occurs between those 
levels as well as between the two levels below n = 6 and 
n = 8. The largest m values are almost unpaired, for 
they would have lead to non-zero total angular momen- 
tum. We conclude that in the low filling regime and at 



intermediate interaction strength we can understand the 
FFLO pairing mechanism in a trapped system as pair- 
ing between fermions from different harmonic levels. We 
observe that the pairs are formed in such a way so that 
the total orbital angular momentum of all pairs is always 
zero, and the orbital angular momentum is minimized 
for each pair. Finally, similarly to the untrapped case 
where the pairs are produced with a finite center of mass 
momentum (vanishing for the balanced case), the FFLO 
state in the harmonic trap corresponds, in a classical pic- 
ture, to pairs whose center of mass is oscillating around 
the minimum of the trap with an amplitude increasing 
with population imbalance. 

Momentum distributions and density profiles at low 
filling of the lattice. Fermion systems with imbalanced 
populations have been realized experimentally in onc- 
and elongated three-dimensional harmonic traps. The 
density profiles of the populations were found to be qual- 
itatively different in the two cases. In three dimensions, 
one observes the formation of concentric shells where, 
for very low polarization, the core is fully paired, i.e. 
zero local magnetization, and the wings are partially po- 
larized [51 [7]- On the other hand, it was observed in 
one-dimensional systems that, for low polarization, the 
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FIG. 12: (Color online). Single particle, (a) and (b), and pair 
(c) Green functions in the harmonic level basis (QMC results) 
for a medium polarization (P=0.22). The total number of 
particles is 26.9. The diagonal pairing has almost completely 
disappeared and the pairing mostly occurs between the levels 
n = 3 and n = 4. More precisely, the strongest pairing occurs 
between re = (4, —4) and re' = (3, 3) and analogously between 
re — (4,4) and k' = (3,-3). There is, in addition, a small 
contribution from the levels re = (4, —2) and re' = (3, 1) and 
re = (4,2) and re' = (3,-1). Note that each pairing corre- 
sponds to a non-vanishing total angular momentum for the 
pair. 



unpolarized fully paired populations are located at the 
edges of the cloud while the core is partially polarized 
[8]. The role of dimensionality in this qualitatively dif- 
ferent behavior has been one of the focus of studies on 
this system. Consequently, the behavior of the system in 
two dimensions is of considerable interest. 

We present here results of our DQMC study of the 
trapped two dimensional system. The presence of the 
trap imposes constraints which make the simulations 
much harder than the uniform case. The number of par- 
ticles should be large enough so that at large P the mi- 
nority population will still be appreciable but not so large 
that the local density in the core regions is close to half 
filling. Another constraint is that the size of the lattice 
be large enough to ensure that particles do not leak out. 
These constraints limit our ability to do simulations for 
system sizes beyond 20 x 20. 

As for the uniform system, the most important indi- 
cator of the presence of the FFLO state is the pair mo- 
mentum distribution. Although the plane wave basis is 
not the natural one in the harmonically confined case, 



FIG. 13: (Color online). Single particle and pair Green func- 
tions in the harmonic level basis (QMC results) for a strong 
polarization (P=0.37). The total number of particles is 27.4. 
The pairing occurs between the levels n — 5 and n = 3 and 
also n — 4 and n = 2, i.e. with total zero orbital angular mo- 
mentum. Still, there is small contribution to pairing between 
K — (5, —5) and k' = (3, 3) and re — (5, 5) and re' = (3, —3). 



0.036 
0.032 
0.028 
0.024 
0.020 
0.016 
0.012 
0.008 
0.004 
0.000 



8 9 10 



FIG. 14: (Color online). Pair Green function in the harmonic 
level basis (MF results) for a polarization P=0.27. The results 
are similar to the QMC results: pairing is maximum among 
the Fermi-levels n — 7 and n = 9 and also among the two 
levels below n = 6 and n = 8. The largest m values are 
almost unpaired, for they would have led to non-zero total 
angular momentum. 



we study the momentum distributions because they are 
of experimental interest. We will show that despite the 
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shortcomings of this language one can still detect the 
FFLO pairing signal this way. In addition, since the 
trap destroys translational invariance it is very useful 
to study the density profiles and local magnetization, 
m(x, y) = pi(x, y) — pi{x, y) where pi (pi) is local density 
of the majority (minority). We start with the unpolarized 





FIG. 15: Momentum distributions of (a) the single particles, 
ni(k x ,ky) — ri2(k x ,ky) and (b) the pairs n pair (k x , k y ). The 
total number of particles is 22.3, P = 0, /3 = 10, U = — 3.5t, 
the trap potential is Vt = 0.065 and the lattice size 20 x 20. 



system. Figure |15[ a) shows the momentum distribution 
of the particles (the two populations are identical) for a 
system with a total of 22.3 particles, P = 0, f3 = 10, 
U = — 3.5t and a lattice size of 20 x 20. The trap po- 
tential is given by V t — 0.065. Figure [l5^b) shows the 
pair momentum distribution and exhibits a sharp peak at 
zero momentum. Now we polarize the system keeping the 
total number of particles constant which corresponds to 
the experimental situation. Figure [i"6] shows the momen- 
tum distributions of the (a) minority and (b) majority 
populations and (c) the pairs. The system has a total 
of 21.4 particles, P = 0.55, /3 = 10, U = -3.5i and a 
trap potential V t = 0.065 on a 20 x 20 lattice. The Fermi 
temperature of the system is Tp — 1.86. Figure [l6[c) is 
qualitatively different from Fig. [l5^b) and shows clearly 
that when the confined system is polarized it exhibits 
FFLO states with pairs forming with nonzero center of 
mass momentum. This behavior was observed for a wide 
range of polarizations and interaction strengths. The ver- 
tical scale in Fig. 16 c) shows that the number of pairs 
is very small. This is due to the small total number of 
particles in the system. A simulation for a larger system 
but with the same characteristic density [37] should give 
a stronger signal in the form of higher peaks at nonzero 
momentum. This effect of the total number of particles 
was shown in the one dimensional uniform case in Ref. [S] . 

During our simulations we measure the density pro- 
files of each species and we calculate the local magneti- 
zation m(x,y) = n\(x,y) — ri2(x,y). The profiles shown 



FIG. 16: The momentum distributions of (a) the minority 
and (b) majority populations and (c) the pairs. The total 
number of particles is 21.4, P = 0.55, j3 = 10, U = — 3.5t on 
a 20 x 20 lattice. The trap potential is Vt = 0.065 and the 
Fermi temperature is Tp = 1.86. 



in Fig. [Tf] correspond to the situation when FFLO-type 
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pairing has been observed in the system as in Fig 
One observes that the system is partially polarized at 
the core and fully polarized in the wings (where we see no 
minority particles). There is no fully paired phase where 
m{x,y) would disappear within the size of the cloud. 

Density profiles are the basic quantities that character- 
ize the trapped system. The first experimental results in 
a three-dimensional system show the formation of con- 
centric shells where for very low polarization the core 
is fully paired (no local magnetization) and the wings 
are partially polarized (see [6] and [7J). On the other 
hand in the one-dimensional system it has been observed 
that there exists a low polarization regime where the un- 
polarized superfluid is located at the edge of the cloud, 
and the core is partially polarized [8]. The issue of this 
dimensionally driven transition caused considerable in- 
terest. It is interesting to look at the intermediate two 
dimensions and study the behavior of the density profiles 
to see whether it follows more closely any of the two lim- 
iting scenarios. During our simulations we measure the 
density profiles of each species and calculate the local 
magnetization m(x,y) = ni(x,y) — ri2(x,y). The pro- 
files shown in Fig. [IT] correspond to the situation where 
FFLO-type pairing has been observed in the system as in 
Fig. 16 One observes that the system is partially polar- 
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FIG. 17: Density distributions of majority (ni(x,y)), mi- 
nority (ri2(x, y)) and the local magnetization (m(x,y)). To- 
tal number of particles is 21.4, P=0.55, /3 — 10, Lattice size 
20x20, Trap potential Vt = 0.065. 



ized at the core and fully polarized in the wings (where 
we see no minority particles). There is no fully paired 
phase where m(x, y) would disappear within the size of 
the cloud. 

In the very low polarization regime we observe oscilla- 
tions appearing in the profile of the local magnetization. 
We looked in detail into these results in order to establish 
whether the oscillations are linked to the FFLO type pair 
density wave behavior. We found however that the oscil- 
lations are present in the system even when there is no 



interaction between particles as seen in Fig. 18 From the 
discussion in the preceding section, where we have shown 
the relevance of the harmonic levels at low fillings, we be- 
lieve that this effect stems from the underlying harmonic 
level structure. In the balanced case, it has already been 
shown that the density of a fermionic cloud in a trap can 
exhibit oscillations with minima or maxima in the cen- 
ter of the trap depending on whether the last filled state 
corresponds to an odd or even harmonic level [39 . 

Harmonically confined system around half filling MF 
study. As mentioned earlier, the Quantum Monte-Carlo 
method suffers from a stronger sign problem for higher 
fillings of the lattice with the trap. However, we success- 
fully studied the system imbalance around half-filling of 
the lattice in the trap using the Mean-Field method. In 
the figures [l9j [20j [2TJ [22] the order parameter is shown 
in real space as well as in Fourier space, for increasing 
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FIG. 18: (Color Online). Cut through the center of the trap 
showing the local magnetization (m(x,y)). Comparison of 
interacting and non interacting profiles for low polarization 
using MF and QMC The oscillations seen in the magnetiza- 
tion are present even in the non-interacting situation (dashed 
line). From both the MF and QMC, one can see that the 
interaction might change the profile, but does not crucially 
change the oscillation pattern. Therefore, we attribute the 
oscillations to the underlying harmonic levels rather than to 
the FFLO order. 



value of the polarization. The numerical results were ob- 
tained for a lattice size 41 x 41, an interaction strength 
U = —5 and chemical potential at the center of the trap 
corresponding to half filling. 

At low polarization (P=0.13), Fig. 
is similar to the balanced case, i.e. 
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a maximum num- 
ber of pairs at the center of the trap, decreasing on the 
border. The Fourier transform simply depicts a peak at 
k = 0, emphasizing BCS-likc pairing. At higher polariza- 
tion P = 0.43, Fig. [20j a structure in the pairing order 
A appears at the center of the trap, leading to clear os- 
cillations in Fourier space. This pattern appears first at 
the center of the trap simply because it corresponds to 
half filling which, as explained in a previous section, is 
strongly unstable towards the FFLO state. Indeed, this 
is emphasized by the two figures [21] and |22| correspond- 
ing respectively to polarization P = 0.48 and P = 0.66. 
The checkerboard pattern of |A| 2 in real space becomes 
more and more visible. Note that similar results have 
been previously shown in |32j . However we would like 
to emphasize the link between this pattern and the na- 
ture of the pairing in the homogenous situation. Indeed, 
in Fourier space four peaks are clearly observed. Their 



positions, (k x = 0, k y = ±q) and (k x = ±q, 



k y — 



0), 



precisely match the ones observed in the homogeneous 
situation, both in the QMC results and in the MF ones, 
around half filling. In addition, one can see that the os- 
cillation period of the order parameter becomes shorter 
with higher polarization, i.e. corresponding to a larger 
center of mass momentum q of the pair, which is de- 
picted by the spreading of the four peaks further away 
from k = 0. This also shows that the oscillations in real 
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space are not related to the underlying harmonic levels, 
but really to the FFLO order. From the experimental 
point of view, this signature of the FFLO order could 
be measured either directly in the density of pairs or in 
their velocity distribution. Of course, the present mean- 
field calculation does not include the thermal fluctuations 
which are crucial to properly describe the condensation 
of the pairs which, at large interaction, arises at a tem- 
perature ksT m t 2 /U lower than the pair formation tem- 
perature k B T K.U [33H36] . 




FIG. 19: Mean Field parameter A as a function of the posi- 
tion (top) and in Fourier space (bottom) for a low polarization 
value (P=0.13), around the half-filling situation, in the pres- 
ence of an harmonic trap. The structure is similar to the 
balanced case, i.e. a maximum number of pairs at the center 
of the trap, decreasing on the border. The Fourier transform 
simply depicts a peak at k = 0, emphasizing a BCS-like pair- 
ing. 



V. CONCLUSIONS 

Our results, based on QMC and MF calculations, 
strongly emphasize that the FFLO state is the ground 
state of the fermionic Hubbard model on the square lat- 
tice for a large range of parameters, both with or without 
harmonic confinement. At low filling, the FFLO state 
is similar to the bulk situation (i.e. particles having 
a quadratic dispersion relation), where the pairs have 
a vanishing total angular momentum, but a finite ra- 
dial component for the center of mass momentum. On 
the contrary, around half-filling, the underlying Fermi 
surface due to the lattice structure, leads to a FFLO 
state having only discrete value of the center of mass 
momentum, namely around (kx — 0,ky — ±q) and 
(kx = ±q,ky = 0). We have given an explanation in 
terms of matching fermionic momentum on the Fermi 




FIG. 20: Mean Field parameter A as a function of the po- 
sition (top) and in Fourier space (bottom) for a polarization 
value P=0.43. A structure in the center of the trap is clearly 
visible, leading to oscillations in the Fourier transform. 




FIG. 21: Mean Field parameter A as a function of the po- 
sition (top) and in Fourier space (bottom) for a polarization 
value P=0.48. The checkerboard pattern is a clear signature 
of the FFLO state. The Fourier transform depicts four peaks 
at the positions (kx — 0, ky — ±g) and (kx = ±q,ky — 0), 
precisely like in the homogeneous situation at half-filling. 



surfaces. We have also shown that, in the presence of an 
harmonic confinement and at low fillings, the harmonic 
level basis gives rise to a simple understanding of the 
pairing mechanism. In addition, we have shown that the 
harmonic levels are at the origin of the oscillations seen 
in the local magnetization, which, therefore, are not a 
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FIG. 22: Mean Field parameter A as a function of the 
position (top) and in Fourier space (bottom) for a polariza- 
tion value P=0.66. The checkerboard pattern depicts now a 
shorter period in real space, translating into a larger spread- 
ing of the four peaks in the Fourier space and corresponding 
to pairs having a larger center of mass momentum compared 
to Fig. [21] 

signature of the FFLO state. Finally, still in the pres- 
ence of an harmonic confinement, but around half filling, 
we have shown that the pairing mechanism is essentially 
identical to the homogeneous situation, leading to clear 
signatures in the pair density, both in real space (checker- 



board pattern) and in Fourier space (four peaks), which 
allows for a possible experimental observation with cold 
atoms. 

In the presence of a harmonic trap, it would be in- 
teresting to study the dynamics of the pairs in response 
to a sudden quench from balanced to imbalanced popu- 
lations where our study indicates that one could expect 
to observe the oscillations of the center of mass in the 
trap. In addition, from a mean field point of view, the 
following points would be interesting to consider. By 
monitoring the wavelength and the amplitude of the os- 
cillations of the order parameter, one should be able to 
determine the nature of the pairing and possible transi- 
tions between paired phases. One should also take into 
account the effects of terms beyond mean field to deter- 
mine properly the critical temperature of the transition 
(BKT-like) and to estimate the strength of the quantum 
fluctuations thus allowing for a better comparison with 
possible experimental results. Finally, one could study 
more exotic situations, like asymmetric tunneling rates, 
or in the presence of an effective gauge field. 
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